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Abstract 



Using Lagrangian formalism we consider evolution of the ensemble of inter- 
im 

acting magnetohydrodynamic cyclones governed with Langevin type equa- 
tions in the rotating medium. This problem is relevant to the planetary 
cores where the Rossby numbers are small and geostrophic balance takes 
place. We show that variations of the heat flux at the outer boundary of the 
^ spherical shell modulates frequency of the reversals of the mean magnetic 



field that is in accordance with the 3D dynamo simulations. Two scenarios 
of reversals were observed. Either the axial dipole decreases in favour of 
quadrupole and then grows in opposite direction or the mean dipole tilts and 
reverses without decrease of its amplitude. 
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1. Introduction 

The present geodynamo models include thermal and compositional con- 
vection which supply the energy to the magn etic field generation in the 



rapidly rotating spherical shells (lOlsonl . 120071 ). Due to the strong non 



linearity of the system and its three-dimensionality the most reliable way 
of study was the numerical simulations. Because of the diffe rence of the 



conve ctive and magnetic typical times in order of magnitudes (IHollerbach 



20031 ) simulation of the geomagnetic field evolution is a very difficult and 
time consuming problem even for the modern supercomputers. The analysis 
of the obtained bulk of data is also very tricky process. This is motivation 
for considering simpler models, which can demonstrate some properties of 
the 3D models and present statistics for a longer time periods where the 
complex models fail. This way helps not only to summarise some already 
known results but also to predict directions of the 3D modeling, which has 
to deal with response of the model on various interconnected parameters. 

The choice of such a model should be based on the modern knowledge 
of the dynamo process in the core and be proved by the observations. The 
specifics of the planetary convection is presence of the geostrophic state, when 
the scales along the axis of rotation are much larger than the scales in the 
perpendicular planes. The corresponding anisotropy also takes place in the 
wav e-space and leads to the inver se cascades even for the pure hydrodynam- 



ics (IHejda and Reshetnyaki . 120091 ). The alongated along the axis of rotation 
primitive convective cells are cyclones (anti-cyclones), in which the flowing 
up (down) liquid rotates in such a way, that helicity of the cell is negative 
in the northern hemisphere and positive in the southern hemisphere. This 
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set of cells produces the mean magnetic field at the surface of the planet, so 
that the cellular structure of the magnetic field is already smoothed and the 
dipole configuration dominates. This v ery s imple and naive approach was 



successfully tested in 



Nakamichi et al 



( 120121 ) where the net magnetic field 



was produced by the ensemble of the primitive spins located at the circle at 
the equatorial plane. The spins rotate in the vertical plane, interact with 
each other and "feel" direction of the planet rotation. The energy is injected 
to the system with the random force which mimics buoyancy sources. Such 
a system produces quite reasonable sequences of the mean magnetic field 
reversals and can be adopted to the regimes with the variable frequency of 
the reversals. Here we want to extend this model, so that the spins could 
rotate in the horizontal plane as well. Transition to rotation of the spins 
in the space makes possible to describe precession of the spins around the 
vertical axis and, as a result, to mimic the magnetic pole wandering between 
the reversals and the fine structure of the reversal itself. In fact, we ob- 
serve two different classes of reversals, when the spins are synchronized in 
the horizontal plane during the reversal or not. The first class correspond 
to the stable dipole field during the reversal and the second to the decaying 
field. The other point is influence of the external fields on the system. Here 
we introduce analysis of the heat-flux modulation of the frequency reversals 
and compare our results with the 3D dynamo simulations. We also consider 
how the core-mantle heterogeneity of the heat-flux can produce the preferred 
meridional band of the magnetic pole migration during the reversal. 
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2. Two-dimensional approach 



We briefly repeat results of 



Nakamichi et al. 



(120121 ) , which is an extension 



of the Ising-Heisenberg XY-models of interacting magnetic spins. For mor e 
details of the history of the problem and classification refer to lStanleyl ( 1l97ll ). 
The main idea of the domino model is to consider a system of N interacting 
spins Si, i = 1...N, in media rotating with angular velocity O = (0, 1) 
in the Cartesian system of coordinates (x, y). The spins are located over 
an equatorial ring, are of unit length and can vary angle 9 from the axis of 
rotation in the range of [0, 2ir] on time t, so that Sj = (sin#j, cos#j). Each 
spin Sj is forced by a random force, effective friction, as well as by the closest 
neighboring s pins Sj_-\ and Sj,+ i. 



Following 



Nakamichi et al. 



(j2012[ ). we introduce kinetic K and potential 



U energies of the system: 



1 - 

*(*) = o£* 



i=i 

N 



N 



(1) 



t/(t)= 7 £(^-S J ) 2 + A^(S i -S i+1 ). 



i=l 



n=l 



The Lagrangian of the system then takes the form £ = K — U. Making the 
transition to the Lagrange equations, adding friction proportional to 9 and 
the random force x, 

d dC dC ■ ex (n . 
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leads to the system of Langevin-type equations: 



9i — 27 cos 9i sin 9i + A cos 9i ^ sin + sin 9i+ij 
sin 9i ( cos 6*j_i + cos 9 i+ i 



T 



(3) 



9 = 9 



N i VN+1 



1 ... N, 



where 7, A, constants. The measure of synchronization of the 

spins along the axis of rotation 



1 - 

M(*) = -X> S W 

i=l 



(4) 



will be considered to be the total axial magnetic moment. 

With the appropriate choice of parameters the simulated sequences of 



M(t) resembles p a 



(jNakamichi et al 



eoma gnetic records of the magnetic dipole evolution 



2012I ). and possess some important properties of the ge- 



omagnetic field: has irregular time intervals between the reversals, exhibits 
short drops and recoveries of the field (the so-called excursions of the mag- 
netic field), the reconstructed 3D magnetic field, based on idea that each spin 
is a magnetic dipole, has similar structure as the observable one. 



3. Three-dimensional spin model 

The natural next step is generalisation of the model to the 3D spins 
rotation. This extension helps to dispose one of the critical disadvantages of 
the two-dimensional model. It is obvious, that nutation displacements in 9- 
direction should cause response of the Coriolis force perpendicular to 9- and 
^-directions. The Coriolis force causes precession of the spins around the axis 
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of rotation O in the horizontal plane, like it happens for the magnetic spins 
in presence of the external magnetic field or for the spinning top. Formally 
this effect leads to the new evolution equation for the azimuthal angle (p. 
Below we consider how we can take into account this effect and compare it 
to the known equations of the rotating bodies and magnetic spins. 

3.1. Precession and Landau- Lifshitz- Gilbert equation 

In three dimensions direction of the spin Sj = (sin 9i cos ip iy sin Q { sin (p i7 cos 0j 
is defined by two angles in the local spherical system of coordinates. As be- 
fore we suppose that origins of the coordinates are located equidistantly at 
the circle of the unity radius at the horizontal plane perpendicular to Q. All 
the systems are obtained by translation in space, so that all corresponding 
axis are parallel. 

First consider approximation when terms with 9, Cp and quadratic in 9 
(p are neglected. Then, using that fact that for the precession-like solu- 
tion the term propo rtional to p cos 9 should be included in the Lagrangian 



flMiltat et all 120021 ). this leads to the i th spin 

Li = <pi cos 9i - 7 (ft • Si) m - li, (5) 
where m is integer. The interaction term can be written as follows: 
Xi = X [(Si ■ S m ) + (Si ■ Si-!)) = 



A 



sin 9i (sin 9 i+1 cos^ - p i+1 ) + sin 9,^ cos(y^ - <^»_i)) + (6) 
cos^j (cos^+i + COS^j-i^ 



Further we consider two cases: i) m — 1 corresponds to ferromagnetics, 
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when direction of the external magnetic field is important and the spins are 
co-directional to fi; ii) for m = 2 the both directions ±0 are equivalent, so 
that the potential energy concerned with rotation is symmetric relative to 
the equatorial plane as it was assumed in 2D case ([3]). As according to the 
paleomagnetic observations there is no preferable polarity of the magnetic 
field and the Lorentz force in MHD equations is quadratic on the magnetic 
field, the latter case is more suitable for our tasks. As we see below the 
first two terms in r.h.s. of (J5J) for m = 1 leads to Landau-Lifshitz-Gilbert 
(LLG) equation, which describes precession of the magnetic spin in the non- 
dissipative medium in the external Zeeman magnetic field equal to ft. The 
third term describes local interactions of the spins. Now we write the La- 
grange equations for 4 independent variables (9, (p, 9, (p): 

d dC _dC dJF dTZi 
Jtd9 l ~d9 l + df i + ~d9~~ ' 

(7) 

d dC dC dT dTZi n 
1 1 = 0, 

dt d0i dipi d<fi d(fi 

where 

Fi=2 yi + sin2 °i Vi 2 ) , TU = — ^ (Pi Xi + <P% ^i) , (8) 

and tp is the random function. 

Substitution of (J5J) and (jSj) in (j7j) for m — 2, leads to the following system 
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of equations: 



6i — k sin 9i (pi 



X' 



sin 6i sin Oi y/r 



-2 7 cos0i + + 



0. 



(9) 



sin 6 1 



sin Oi sin ^ -^/r 



The exact form of the terms with derivatives of Zj are given by equations: 



X' 



sin Oi 



sin^+i sin(v?i - (p i+1 ) + sin0i_i sin(y?i - ^ 



i-lj 



X' 
sin Oi 



Like before in 2D case ([3]), we use the periodical boundary conditions: 



cot Oi y sin i+1 cos(<p» - ip i+1 ) + sin cos(<£>» - </?i-i) J 

COS^j+i + COS#i_i 



(10) 



O = N , N+1 = 1 , if = tp N , tp N+ i = (ft, i = 1 . . . N. (11) 

Since Lagrangian (jSJ) does not include quadratic terms proportional to 2 , 
ip 2 , equations fl9]) present balance of the forces. Before we come to analysis of 
the system (|9]), we consider case m = 1, where the term with 27 cosO in 
changes to 7. Puting A = 0, k = 0, e = 0, yields = 0, <p = 7, i.e. precession 
around axis z, that corresponds to the solution of the LLG equation without 
dissipation for ferromagnetics: 



Si = - 7 s, x n, 



(12) 
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which has integral of motion — = 0, and for the constant in time f2 the 
d 

other integral: — (Sj • f2) = 0, that follows to 9 = 0. The x-component of 
equation (fl2|) gives the already mentioned above equation for the azimuthal 
angle: ip = 7. Existence of precession distinguishes the 3D case from 2D, 
where the role of rotation comes to attraction of the spins to the poles. ([3]), 
so that derivative moves from ^-component equation to (^-equation. 

For the case m = 2 one gets precession-like equation 9 = 0, up = 2 7 cos 9, 
which predicts reversal of the angular velocity of the spins at the equator 
plane. Obvious, that due to the distinguished direction concerned with ro- 
tation the whole system is not symmetric to the change of z — > —z and the 
break of the reflection symmetry for m = 1 gives preferable polarity of the 
magnetic field (sin 9 7^ 0). For m = 2 one has sin# — > 0, however (p changes 
its sign at the equatorial plane z — 0. Below we consider this phenomenon 
in more details. 

It is instructive to consider balance of the curvilinear terms and the Cori- 
olis force in Vg- and V^-components of the Navier-Stokes equation in the 
spherical system of coordinates with V r = 0: 

-I/ 2 co%9 = H V v cos 9, V e V v cot 9 = -HV e cos 9, (13) 

where H is the amplitude of the Coriolis force, and tangential velocity (V$, V v ) 
is (9, ip sin#). The both equations in (fl3|) gives ip = —H, that corresponds 
to the case with m — 1 and Vg remains undefined. As we see, condition 
of equiprobability of the direct and inverse polarity for m = 2 changes the 
whole scenario. 
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Return to (JH]) with m = 2. Since we have we omitted quadratic in veloc- 
ities terms, this approximation is valid for the slow regimes with precession 
between the reversals of the field. Anyway we consider how the predicted 
asymmetry of the precession velocity is influenced by the random force. Even 
for the small random force distribution of V V (M) in FigJTh^ for the regime 
without reversals in Fig{Tb has zero mean value. The reason of such a seem- 
ing discrepancy with our estimate dp ~ cos 9 is the following. Consider aver- 
aged in time coupled equations in the simplified form without interaction, 
diffusion and random forces for small 0-^1 and constant in time between 
the reversals of the magnetic field 6 = and check if do still changes sign 
for the coupled system: 9 ~ 6 do = 0, 6 do + 6 ~ 7(1 — 6 2 )6 that leads to the 
contradiction: 9do ~ 76 1 7^ 0, that means that the other terms should be 
included in the consideration and naive prediction of change of the sign for do 
is wrong. We conclude, that the random forcing makes system symmetric to 
the change of (^-direction. Here we do not consider extremely small e regimes 
to preserve the realistic for the geomagnetic field level of the field fluctuations 
between the reversals. 

The considered approximation with the small fluctuations of 9 is reliable 
for some astronomical problems without nutation, as well as to the ferro- 
magnetic problems when the temperature is less than the Curie point and 
^-variations are small. This review was used from the methodological point 
of view and now we come to the full equations applicable to the highly non- 
linear regimes with the large accelerations. 
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3.2. Precession and nutations 

Now we consider approximation with large 9 . Cp\ 9 2 , cp 2 using ana. 



ogy to 



2005) 



the spinning top with unity moments of inertia (ILandau and Lifshitz 
with Lagrangian in the form: 

d=~ (e 2 + sin 2 Bi (p\ + {(fi cos 6i + Ci) 2 ) - 7 (O • S,) 2 - X, + * i? (14) 

where Q is angular velocity of the top. In our case ( is fixed and given 
property of the spin. For simplicity we take it equal to unity: C 2 = |Sj| 2 = 1. 
The new variable \I/ is the potential used latter for description of the external 
forces. 

Finaly Lagrangian can be written as follows 
1 • 1 

C i = o e i + 2 $ + ^ cos ^ ~ 7 cos2 9i ~ Ii + **» ( 15 ) 

where constant is omitted because Lagrange equations includes only deriva- 
tives of £. Neglecting of the quadratic terms in velocities in ( 115)) leads to the 
simplified version of Lagrangian (JSJ). 

Substitution of (|T4|) in ((7j) yields to the dynamic equations: 

0< + ifi sinOi - 7 sin 20; +l' ie + k9\ + ^L- V' ig = 0, 

eh ^ 
Cp - sin 9i 9i +X[ + k sin 2 9i <pi + - ^' i(p = 0. 

Neglecting the second in time derivatives and quadratic terms leads to the 
original system (J5J). 

Note, that in contrast to 2D case we already have two equations, that 
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requires "energy" balance fulfilment, and the new numerical methods should 
be used. This point is enforced by that fact, that there are terms with sin8 
in the denominator. That was the reason to use Newton-Raphson iterative 
method in the residual form described in the Appendix. 

Now we consider important for the geomagnetism question on the pre- 
dictability o f the geomagn etic field reversal on the secular variations before 



the reversal (Jacobs 



20051 ) and take as a measure of variations the amplitude 



of the magnetic dipole precession velocity V v around the geographical pole. 
Using equations (JIB]) we generate 30 reversals of the field, see Figf5]^l). Taking 
the typical time t{ between the reversals as 3 10 5 y we come to the estimate of 
the full time interval as 9 10 6 y and the time unity r u = 3 10 5 /2 10 4 = 4 500y. 
The estimate of the azimuthal velocity V v = 0.2 leads in the order of mag- 
nitude to the typical archeomagnetic time estimate of the pole wandering 

7T 1 

(similar to the west drift velocity) r a = r u ~ 3 500y, where we took into 

8 V v 

7T 

account, that between the reversals the pole locates in the cone 6 < — . 

P 8 

We inspect the set of the five reversals of the magnetic field, see Figj2]^2). 
The particular feature of this regime is existence of the intermediate state 
corresponding to the small \M\ when the system does not know in what 
direction it should go, and as a result some of the reversals fails producing 
excursions, see Figj2]^l). To check possibility of the reversal predictability we 
consider the mean values of M and \V V \ over these five reversals, see Figj2|3, 
4). As we see, the change of polarity does not lead to the substantial change 
of |V^|, that means that the behaviour of the model is symmetric in respect 
to the moment of the reversal. 
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4. Magnetic field of the spins 

The essential disadvantage of the domino model based on the assumption 
of the long-range interactions for the spins is absence of connection between 
the spin's position and its magnetic field. That was the reason why we still 
used integral characteristic M as a measure of the field. Let us summarise 
requirements to the model which we want to satisfy: i) the potential energy 
of the spins should have minimum when they are located near the poles and 
the spins should have the same direction; ii) the model provides calculation 
of 3D magnetic field at least at some distanc e from the spins. I dentifi cation 
of the spins with the magnetic dipoles helped 



Nakamichi et al. 



(12012( 1 to get 



3D distribution of the magnetic field for the spins in the form: 

B= 3r(d-^-r»d (17) 

where d is a magnetic moment of the dipole in the centre of the coordinates 
and r is the radius-vector of observation. It should be noted, that as follows 
from (IPTjh the stable state of two magnetic dipoles at the distance R with 
corresponding potential energy of interaction 

3(d 1 -R)(d 2 -R)-(d 1 -d 2 )i? 2 , , 

u = ^ ( 18 ) 

is quadrupole, when the magnetic moments di, d 2 are anti-parallel. In its 
turn this would correspond to the death of the dipole magnetic field at large 
R and relation (TPS1) is out of interest for geomagnetism. 

To overcome t his pr oblem and to get the self-organised mean dipole field 



Nakamichi et al. 



( 120121 ) suggested to use the simplified version of (fT51) with 
13 



omitted first term and opposite sign for the second term, see ([T]). This form is 
similar to that one in the ferromagnetics where magnetic domains enforce the 
external magnetic field choosing the same with it direction. In ferromagnetics 
this effect has the pure quantum origin c oncerned with the exch ange energy 



of the magnetic domains. As regards to (iNakamichi et all 120121 ) it was just 
postulated as a convenient way to get similar to observations results and we 
join this idea in the rest of the paper. 

Using equation ( TT7|) for the set of N magnetic dipoles, located at the circle 
of unity radius we compute the vector magnetic field, see distribution of its 
radial component at the distance of the three unity length TZ during one of 
the reversal in Figj3j One of the specific features of this reversal is existence 
of the preferred meridional band where the reversal occures, see FigJ3f2). 
To consider this effect in more details we come to analysis of the spectral 
properties, see evolution of the axial dipole Gauss coefficient g®, in FigJU^l), 
at the distance TZ for the five reversals in FigJ2^1). We conclude that the 
typical time of the drop and recover of the field during the reversal in this 
model is the same, there is some increases of the field just before and after 
the reversal. So as the amplitude of the spins is fixed the net magnetic flux 
of |B| is constant in time as well as its spectrum, during, e.g., the reversal. 
This leads to the exchange of the magnetic energy between the dipole and 
quadrupole modes in FigIJj2), where the ratio of the dipole and quadrupole 
modes gives T>/Q « 13 compared in order of magnitudes with that at the 
surface of the Earth. In the moment of the reversal t ~ 2 100 T> drops, and 
intensity of the qudrupole Q increases to the level of the dipole field. This 
effect is known in 3D simulations as well. The reason of such a distribution 
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is that magnetic energy in the geostrophic systems is quite large and hardly 
can be immediately transferred to the kinetic energy during the reversal. It 
means that the drop of the dipole field should be compensated with increase 
of the higher modes. 

The other, may be the more specific feature of the domino model is that 
for the other four reversals the amplitude of the total dipole field during 
reversal is constant. In contrast to the previous case, where due to the 
uncorrelated horizontal projections of the dipole field the averaged dipole 
field decreased, the second case corresponds to the coherent state when even 
during the reversal all the spins have the same direction. For this scenario 
the mean dipole field rotates in the meridional plane without decrease of the 
amplitude. As we see, this three-dimensional phenomenon for the considered 
regime is more expected and caused by interaction of the spins (fT6|) in in- 
direction during the reversal with the short interaction time compared to the 
times of the reversal and precession. The both scenarios does not contradict 
to observations. 

5. Thermal flux heterogeneities 

One of the important results of the geodynamo theory is that frequency 
of the reversals depend on the spatial d i stribu tion of the heat flux at the 



core-mantle boundary (IGlatzmaier et al. 



19991 ) . Fluctuations of the heat 



flux of order 10-20% of the mean flux value caused with the processes in 
the mantle and D" layer have typical time scales 10 6 -10 7 y are much larger 
time scales in the liquid core 10 4 -10 5 y. In particular increase of the heat 
flux towards the poles is equivalent to increase of rotation, that damps the 
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reversals. Such fluctuations appears to be the thermal traps for reversals. In 
its turn, the decrease of the heat flux in the high latitudes leads to chaotic 
behaviour of the magnetic dipole and increase of the reversal frequency, when 
the Archimedean forces become more significant. This effect corresponds to 
increase of the Rossby number, which has threshold at 0.12 where regime 
of the stable dipole field changes to the regime of the frequent reversals. It 
looks tempting to reproduce this effect in the domino model, where the large 
number of the reversals can be easily simulated. 

We now extend the concept of the spin from the purely magnetic system 
to the whole cyclone system, including its hydrodynamics, and we introduce 
correction <p) to the potential energy U, which takes into account the 
heterogeneity of the thermal flux. The new effective force Fj is proportional 
to the corresponding derivatives of with the opposite sign will appear in 
(|T6]) . Here we study influence of the various forms of \I/ on the behaviour of 
M(t). 

Let 9, ip) = C^tpit, 9, <f), where C$ is a constant, and the spatial 
distribution of the potential is given by ip = — cos 2 9. Accordingly to the re- 
cent estimates of th e heat-flux variatio ns at the core-mantle boundary, which 
can be about 20% ( Olson et al. . 2002 ). we conclude that Cw, ~ 0.2-^= ~ 1. 
Then, > corresponds to the stable state in the polar regions, 9 = 0, tt, 
and the appearing force F = — sin26 l , acting on the cyclones, is directed 
towards the poles. This regime corresponds to the increase of the thermal 
flux near the poles that causes the stretching of the cyclone along the axis of 
rotation. So as F changes sign at the pole 9 = tt, and the force is directed 
to the pole (or outward) for 9 — > tt — and 9 — > tt + 0, we can consider it 
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as the spherical coordinate when it is needed. In FigJS]we demonstrate the 
effective influence of F on M for the regime in[2](l). The increase of the 
thermal flux along the axis of rotation leads to the partial suppression of the 
reversals of the field, Figj5^1). Note that, for the chosen potential barrier ip, 
the dependence of F is equal to the 7-term in ([3]): the increase of the thermal 
flux at the poles leads to the effective increase of rotation and amplification 



of geostrophy, caused by the rapid daily rotation of the p 



are in agreement with the 3D simulations, see Fig. Id in (IGlatzmaier et al. 



anet. Our results 



19991 ) . The further increase of C$ = 2) leads to the total stop of the 
reversals. It is more interesting that, using even larger > 10, one arrives 
at regimes with a nearly constant in time \M\ < 1 defined by the initial dis- 
tribution of Sj. In other words, the super flux at the poles can fix th e spins 



which are still not coherent. There is some evidence (jShatsillo et al.l . 120051 ) 
that the geomagnetic dipole in the past could have migrated from the usual 
position near the geographic poles to some stable state in the middle lati- 
tudes. Within the framework of our model, we can explain this phenomenon 
by the thermal super flux at the poles. Later we will discuss some other 
scenarios which yield similar results. 

For negative C^, when the geostrophy breaks due to the relative in- 
tensification of convection in the equatorial plane, we get the opposite re- 



sul t, see FizM2): t 



r e reg ime of the frequent reversals observed in Fig.lc 



in (IGlatzmaier et al. 



1999j). In this case force F is directed from the poles 



and the equilibrium point at the poles becomes unstable. The new mini- 
mum of the potential energy at the equator leads to the appearance of a 
new attractor, with the dipole at the low latitudes. Similar behaviour of the 
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mag netic dipole is obse rved on Neptune and Uranus; for more details see, 



e.g. (jCupal et al. 



2002|). 

Now we consider the non-axial-symmetrical potentials. Let ip = cos 2 9 ■ 
cos 2 ip. It is clear to see, that for the large |C^| frequency of the reversals 
should follow the previous axial-symmetrical case of ip, in particular for = 
3 the number of reversals decreases, see Figj5](3). Moreover, this form of ip 
with dependence on <p leads to the new effect: appearing of the preferred 
meridional band of the magnetic poles migration during the reversals. We 
compare distribution of h\{g\) without potential, see FigJSfl), with that one 
with the non-axialsymmetrical potential. For the first one we have isotropic 
distribution of the horizontal dipoles, see Figj6^1). Insertion of the potential 
leads to appearing of the barriers for the reversals with gaps at ip = and 
(p = 7r, see Figj6](2), where the large circles corresponding to the magnetic 
dipole at the equator plane have the elliptical distribution, and near the poles 
(small circles) distributed isotropicaly. For the considered form of ip sign of 
Fg did not change in (p. 

The last example is with F$((p) changing sign at the equator plane: 
if> = cos 2 6 ■ sin2<^, see Figj6]^3), with the only one reversals for the model 
simulation and the anisotropy observed for all positions of the integral dipole, 
even for the locations near the poles, when the amplitude of horizontal dipole 



\f g\ 2 + h\ 2 is small. In other words, fluctuations of the heat flux causes not 
only preferred meridional bands of the magnetic poles migration, but the 
anisotropy of the poles locations between the reversals as well. We also ob- 
serve states, when for the large |C^| reversals stop because some of the spins 
come to the thermal trap near the poles (for some particular values of ip) and 
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the full reversal which includes reversals of all spins is impossible. 
6. Conclusion 

It appears, that even a small extension of the original toy domino model 
leads to some nice physical effects known in the various fields of physics. In 
spite of simplicity of this approach it has very strong advantage: one can test 
very different scenarios of magnetic field evolution basing on some realistic 
information on the cyclone convection in the core. This approach is well 
designed for modelling of interaction on the different scales and number of 
the neighbour cyclones, which can be obtained from the 3D models. The 
other possibility is to use two distinct sets of cyclones inside and outside of 
the Taylor cylinder, which separates two different regions in the liquid core. 

We also believe, that domino models can be modified for the analysis 
of the net kinetic and may be current, magnetic helicities as well. The first 
kind of helicities is the most preferable for the analysis, because the cyclone's 
hydrodynamics to some level is stable during the reversals, what is not the 
case for the generated in the cyclone magnetic field. The latter obstacle 
makes interpretation of the results in the domino model not so unambiguous. 

A. 

We rewrite system (Tl6|) in the form of the first order differential equations: 

Vi-6i = 0, 
Wi-<j>i = 0, 

V + kV + Wsm0i + A = 0, 
W - sinflj V + k sin 2 6 t W + B = 0, 
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(19) 



relative to the vector 



where 



W n 

\r J 



(20) 



A = - 7 sin2^+j; e + i|-^, B =%+ e -j=--^6^- (21) 



Using the implicit first order Euler scheme in time for the n th time step we 
introduce residual vectors for ffT9l): 



( V n dt - e n + Q n ~ x ^ 
W n dt-ip n + ip™- 1 

yn _ yn-1 + K yn df + yyn ^ Q n ^ + ^ 
W n _ W n-l _ gin Q n yn ^ + R ^2 g B jyn ^ + fi n ^ 



(22) 



The iterative Newton-Raphson process for the p -iteration gives: 



yp y^p— i 



(23) 
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where Jacobian matrix has form: 



J 



<9e™ 



dt 


1 + ndt 





dt 

sin 6? dt 
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sin 0? <ft 1 + « sin 2 0? dt 



1 3 



-1 



cos 0?dt+^dt ^-d/ 



dtp n 



dt 



(24) 



J 4 3 = _ cos 0" y n dt + k sin 20" W n dt + 



35™ 
30^ 



dt. 



(25) 



For each time moment n and spin z one has iterative process ff23l) with up- 
dated other spin values. After the desired convergence is reached transition 
to the new time step n + 1 is done. At the each time step we check conditions 
0j G (0, 2 it), (fi G (0, 2 ir) and correct values if necessary. This algorithm 
appears to be stable to singularity at the axis and can be easily transformed 
to the second-order Crank-Nicolson scheme in time. Anyway this problem 
is more stable than the PDE in the spherical coordinates, because the large 
derivatives near the pole forces spins to go out from the pole region and 
instability in such a self-organised way is damped. The spins are squeezed 
from the pole regions. 
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Figure 1: Dependence of velocity precession V v on dipole polarity M (1) and evolution 
of the dipole M in time (2) for 7 = -1.8, A = -6, k = 0.1, e = 0.03, r = 0.01. 




Figure 2: Evolution in time of M for 7 = 1, A = —4.8, e = 0.8, /c = 0.2 (1) and reversals 
of the magnetic field (on degree of thickness of the line) for 5 time intervals (9640, 9690), 
(12000, 12050), (12260, 12310), (17840, 17890), (18670, 18720) where for the 3rd and 5th 
interval sign of M is changed (2); (3) - the mean value M for the previous plot (2) and 
the mean precession velocity | V v | vice M (4) . 



Figure 3: The map of _B r -component of the magnetic field in the Mollweide projection for 
the time moments t = 2500 (1), 2600 (2), 2700 (3). The white colour corresponds to the 
positive values and the black one to the negative. 




Figure 5: Evolution of M for the regime in Figl2ji with ip = cos 2 6: = 2 (1), CV, = — 1 
(2) and — 3 with ip — cos 2 9 cos 2 <p (3). 
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